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Hot, Rotating Disks In General Relativity: Collisionless Equilibrium Models 

A. Katrin Schenk , Stuart L. Shapiro 2,3 , and Saul A. Teukolsky 1,4 

ABSTRACT 

We present a method for constructing equilibrium disks with net angular momentum 
in general relativity. The method solves the relativistic Vlasov equation coupled to 
Einstein's equations for the gravitational field. We apply the method to construct disks 
that are relativistic versions of Newtonian Kalnajs disks. In Newtonian gravity these 
disks are analytic, and are stable against ring formation for certain ranges of their 
velocity dispersion. We investigate the existence of fully general relativistic equilibrium 
sequences for differing values of the velocity dispersion. These models are the first 
rotating, relativistic disk solutions of the collisionless Boltzman equation. 



1. INTRODUCTION 

Rotating stellar disks in dynamical equilibrium have a long history in astrophysics. Such 
systems are described by self-consistent solutions to the Vlasov equation for the phase-space 
distribution function / coupled to the equations for the gravitational field. Even in Newtonian 
gravitation, finding solutions is difficult since the configurations are nonspherical and have a 
relatively large number of nontrivial phase space degrees of freedom. (For a general review and 
discussion, see Fridman & Polyachenko 1984 or Binney &: Tremaine 1987| .) When the configuration 



is relativistic, the gravitational field is described by Einstein's equations and hence the problem is 
even more difficult. 



In a previous paper (Shapiro & Teukolsky 1993; hereafter [Paper I|) , we developed a method 
for constructing equilibrium axisymmetric star clusters with net rotation in general relativity. 
In this paper we adapt this method to treat disks. Once again we restrict our attention to the 
simplest phase-space distribution functions that can give rise to nonspherical equilibria, functions 
of particle energy E and angular momentum J z alone. Because E and J z are integrals of the 
motion, choosing a distribution function of the form / = f(E,J z ) guarantees that we will have a 
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solution to the Vlasov equation, provided the metric is determined self-consistently. No further 
dynamical equations need to be solved for the matter. By contrast, in equilibrium fluid systems 
one must integrate the equation of hydrostatic equilibrium. 

Several researchers have studied the properties of self-gravitating fluid disks in general 
relativity. Morgan &: Morgan (1969)| presented an analytic model of a fluid disk with no net 



angular momentum. This work was followed by Bardeen &; Wagoner (1969) who studied rapidly 



rotating, semi-analytic, zero pressure models and by [Salpeter fc Wagoner (1971)| who looked at 



rapidly rotating fluid disks with some thickness. Recently the model presented by Bardeen & 



Wagoner (1969) was solved analytically by Meinel (1997)| . Collisionless disks with no net angular 



momentum have been modeled and evolved using a particle simulation code by Abrahams et.a\ 



(1994) . The work presented here presents the first fully general relativistic, rotating, collisionless 



disk models. 

As an illustration of our method, we will focus on the relativistic generalization of an 
important class of Newtonian disks, the Kalnajs disks ( [Kalnajs 1972 ; Binney & Tremaine 1987). 



These disks are completely described by simple analytic expressions. In the limiting case that 
all the particles move in circular orbits, the angular velocity and surface density are just those 
of the corresponding fluid Maclaurin spheroid in the disk limit (eccentricity — * 1). In Newtonian 
theory, however, it is known that equilibrium disks supported against collapse by rotation alone 



are unstable to ring formation (see, e.g., Binney fc Tremaine 1987 , §5.3). Kalnajs (1972) showed 



that the disk can be stabilized by "heating" it, that is, converting some of the ordered rotational 
energy into random thermal motion while keeping the surface density the same. (Properties of 
Newtonian Kalnajs disks are reviewed in Appendix A.) We construct relativistic generalizations 
of Kalnajs disks and study their properties. 

In addition to its potential astrophysical significance for generating models of highly 
relativistic stellar disks or collisionless particle distributions, the method presented here provides 
a class of rotating equilibria that has not been treated previously in the general relativistic 
literature. These collisionless models join Kerr black holes and rotating fluid stars and fluid 
disks as physically realistic rotating equilibria in general relativity. Such solutions to Einstein's 
equations provide important insight into the effects of rotation in a strong gravitational field. 
They also provide useful initial data for evolution codes in general relativity. In particular, such 
codes can assess the stability of relativistic Kalnajs disks, since, as we will see below, we can only 
make heuristic statements about their stability from the results of this paper. 



2. Basic Equations 

In this paper we consider rotating equilibrium stellar disks that are axisymmetric. The metric 
can then be written in the form 

ds 2 = -e^ +p dt 2 + e 2a (dr 2 + r 2 d9 2 ) + e^~ p r 2 sin 2 9(d<f> - udtf (1) 
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where the metric potentials p, 7, u, and a are functions of r and 6 only. Here and throughout we 
set G = c = 1. This is the same form of the metric used in Paper I. 

All calculations in this paper are performed in the ZAMO frame (Bardcen 197C; Lightmar] 
et al. 1975| ). This is the orthonormal frame of a zero angular momentum observer, whose basis 



1-forms u) a are related to the coordinate basis dx a by 

lo & = L & a dx c 



where 



e" 







e v 
re h 



-tor sin 9e@ 







r sin Be 13 



(2) 



(3) 



and where 



P 



1 + P 
2 ' 
7- p 



The corresponding orthonormal basis vectors are given by 





2 r«. 

c a — u a 



dx a 



(4) 
(5) 

(6) 



where L a a is the inverse matrix of L c 



Following Komatsu et al. (1989) , |Cook et al. (1992] , and Paper I, we write the Einstein field 
equations that determine p, 7 and to in the form 



V 2 [pe^ 2 ] = S p (r,p), 



V 2 + -d r 

r 



dJ[^ 2 ]=S,(r,p), 



(7) 
(8) 
0) 



where V 2 is the flat-space, spherical coordinate scalar Laplacian, p = cos 6, and S p , and S u are 
effective source terms that include the nonlinear and matter terms. In the equations below, we 
explicitly exhibit the disk nature of the matter by defining a disk stress-stress energy tensor, i^(r) 
such that the full stress-energy tensor, T^(r,p) is given by 



4(r,p)=4(r)^±. 



(10) 



Thus the integrals used to calculate the stress energy components, equations ([3l|)-(|3"3]) below, 
need only be evaluated in the disk plane (p = 0). 
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The effective source terms are given by 

S p (r,fi) = Rp{r,{i) +t p (r) 

Sj(r,fi) = Rj(r,fi) +ty(r) 
Su>(r, A«) = Ru(r, At) + ^( r ) 



r 
r 



where we have defined matter dependent and matter independent source terms such that 
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=7/2 Jl _ii 
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+r 2 (l- M 2 )e- 2 Ma; 2 r + 



1- M 
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Rj(r,n) 



87re 7/2 e 



2 " 4 (*?+*!) 



=7/2 7 



_1 2 _ !~At 2 2 
2 7 ^ 2r 2 



= 87re^ 2 e 2CT (l + J) (4 + ^), 



e (7"2p)/2 w 



~\ ( 2 ^V + ^7,r) + £ (2p,„ + ^) + ^P 2 r - 7 2 ) 



- T 2 .) - r 2 (l - At 2 )e" 2 ' U + i-/^ 



-<-'l)+|«+0- r(1 _, 2)1/2 



(11) 
(12) 
(13) 



(14) 
(15) 
(16) 
(17) 

(18) 
(19) 



Here tp are the orthonormal components of our disk stress-energy tensor for collisionless matter in 
the ZAMO frame (see below). 

The fourth field equation determines a and is given by 



2-) -1 



- 2 (AV + 7m) " {(1 " At')(l + n, r y + [At - (1 - Ull A ) 
\^[r 2 {l,rr + 7 2 ) - (1 - At 2 )(7, w + 7 2 M )][-At + (1 - ^hj 



+n,r 



1 1/ 2x 

-/i + /w7, r + -(l-/z )7 )A1 



+ 27/J-At 2 + At(l " At 2 )7,^] 



-r(l - At 2 )(7,r M + 7,r7,^)(l + r7, r ) - ^r 2 (p, r + 7,r) 2 
-^(1 - At 2 )(P,r + 7,r)(AV + 7,m) + 4MI - At 2 )(p^ + 7,^) 2 
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— (1 - H )% r {p,r + %r)(ft tl t + 7,/i) 
+ - ^h> 2 (/V + 7,r) 2 " (1 " M 2 )W + 7„) 2 ] 



+(l-/x 2 )e 



-2p 



4 2 i 3 / I 2\ 1 2 /i 2\ 2 

-r /iw r + -r (1 - fi )w >r w iA4 - -r /i(l - ^ )w 



+ ir 4 (l - M 2 )7,r^,r^^ - jr 2 (l - fi 2 hJrW r - (1 - /i 2 )^] 



-r a ^-(l-M 2 )7 rf Je aff M<?-<S)- W 



r 



+r 2 (l - / , 2 ) 1 /2(l + r7jr ) e 2a 87rt rM|_ (2Q) 

The terms containing A and t r ~ in the above equation are identically zero for an equatorial 
disk since each particle has no ^-component of momentum. The remaining matter term in 
equation (p0|), 

-r 2 ^-(l-^ 2 ) 7 ,,]e 2 -4<^ (21) 



r 

also does not contribute since 7,„ L = q vanishes by symmetry. 

The stress-energy tensor for the matter is determined by the phase-space distribution function 
/, which is governed by the relativistic Vlasov equation (the matter dynamical equation). Any 
distribution function of the form / = f(E, J z ) is an equilibrium solution of the Vlasov equation in 
axisymmetry. Here E and J z are the two constants of motion associated with the Killing vectors 
d/dt and d/dcp. For our two-dimensional disk system they are 

d 

E = -p- — = e'V + weT'i-p*, (22) 
at 

J s = p~ = eP'rp^, (23) 

where p a are the orthonormal components of the particle 4-momentum p in the ZAMO frame and 
we have defined, for convenience, the quantities, 

u e = u(r,0), 
(3 e = /9(r,0), 

uj s = w(r, 0), etc. (24) 

Physically, E is the conserved energy of a particle and J z is the conserved angular momentum 
about the symmetry axis. 

Working in two-dimensions, we define the disk stress-energy tensor for the matter, 

t& P = f tfjfii*. (25) 

J p l 
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Here, for a disk, we have 

pi= [( p r)2 + ( p 0)2 +m 2]l/2 j (26) 

and 

d 2 p = dy f dp\ (27) 

where m is the particle mass. For the field equations ([?]) - (|9|) and (pop, we need only the 
combinations 

t X = t\- if, (28) 
i 2 = 4, (29) 
t 03 = 4- (3°) 

Using equations (p5|) - (CT), we get 



' ' [(/)2 + (p^)2 +m 2]l/2^ 

t 2 (r) = / / dp* - ^7rJ( E i (32) 



[(/) 2 + (p*) 2 +m 2 ]V2' 
*os(r)= / d/ / dp^f(E,J z ). (33) 



At each spatial point (r,fj,), the limits of integration in equations ( |3l| ) - (|33|) will be determined by 
the distribution function, as discussed in the next subsection. 



2.1. Distribution Function 

The Newtonian Kalnajs distribution function is given by equation ( |A4[ ). We construct a 
relativistic generalization by letting En — > E — m: 



K 



-1/2 



H>o (34) 



2{(bVJ z /R m ) -(E- m)} - 1/ 2 (1 + b 2 ) 
[-..]<0 

Here V, b and R m are constants which in the Newtonian limit have the following interpretations: 

_ Mean angular rotation rate _ Q QR 



Angular speed of a circular orbit £l c irc V 
and Rm is the matter radius. Note that in the Newtonian limit we have 



(35) 



E-m^E N Jz^JzN, (36) 
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where En and J z at are the Newtonian energy and angular momentum as defined by equations (A5) 
and QA6| ). Thus we recover the Newtonian distribution function /jv, equation ( |A4| ), from 
equation (|34"|) in the Newtonian limit. 



Since the distribution function is nonzero only for positive values of the argument of the 
square root, we can solve for the limits of integration for the stress-energy integrals (eqs. (^Tj) 



-(|33D) by using the expressions ((22;) and fl23| ) for E and J z in equation (34) and setting the 
argument of the square root to zero. Thus we find the limits to be 



and 
with 

and 



T 

Pmax 



(pV) 



0< / <ptnax{P*,r) 

P—{r) < < p+(r) 
a 2 + d 2 -\ 



P± 



l-d 2 
ad 
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1 1/2 



(a 2 + d 2 -l) 1 / 2 



l-d 2 l-d 2 
Here a and d are functions of r only and are given by 



a(r) 



-u e (r) 



1 



and 



d(r) 



e j3 e (r)-v„(r) r 



2 v 



V 
>-5- 



uJAr 



(37) 
(38) 

(39) 
(40) 

(41) 
(42) 



Note that in Paper I the distribution functions considered contained a parameter E max that 
was less than m for bound systems of finite extent. The limits of integration were determined 
by E < E max . Here, one can verify that the limits (|37|)-(|38|) never violate the condition E < m. 
Curiously, equality can be attained even in the Newtonian limit: En — > at the outer edge of a 
"cold" Kalnajs disk. 

The total mass-energy of the system is 



M = -J (2T^-d^T)^ t) d% 



J (-2TI + T 



-gd 3 x. 



(43) 



Here, £^ = d/dt is the time Killing vector. Transforming the integrand to the ZAMO frame, and 
integrating out the <5-function in fj,, we obtain 

M = 2vr y [ti(r) + t 2 {r) + 2uj e e- pe rt 03 {r)] e 2cTe+ ^rdr. (44) 

As we will see in Section ||, this equation plays a crucial role in the iterative solution of the 
combined matter and field equations. 
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3. Diagnostic Probes 

There are a number of useful quantities that characterize an equilibrium system once a 
solution has been obtained. The total angular momentum is given by 

j = |r^)d% 

= / T^V=9d 3 x. (45) 
Transforming to the ZAMO frame, and again integrating out the delta function in ju, we get 

J = 2tt J t 3(r)e 2(CTe+/3e) r 2 (ir. (46) 
The surface rest mass density is given by 

*a(r) =mjdp f J dp^f GR (E, J z ), (47) 

and the total rest mass by 

M = 2ir J a (r)e 2ae+l3e rdr (48) 
The binding energy of the system is defined as 

E b = M - M, (49) 
and can be computed from equations (44) and fl48|). 



4. Scaling and Non-dimensional Units 

The quantities m and M can be scaled out of all the above equations. For example, we can 
define 

p =p/ m , f = r/M, J=J/M 2 

t aP =t a/3 /M~\ f GR = f GR /M- l m-\ J z = J z /mM. (50) 

With these definitions, all the previous equations can be written in tilde variables without m or M 
appearing. Equivalently, the original equations can be solved setting m = M = 1 and scaling the 
final results according to equation ( |50|) to accommodate arbitrary values of m and M. Henceforth 
we will make this simplification. 

In the Newtonian limit there is an additional scale freedom in that we can also set R m /M to 
one. This is not the case for relativistic systems. 
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5. Numerical Solution 

The numerical scheme used here is a 2-dimensional analogue of the procedure adopted in 
Paper | . In this procedure we start with initial guesses for the metric potentials, p, 7, u, and a (see 



section y below). We then integrate equations (|3l|)-(|33|) to find t\, t 2 and io3 up to the constant 
factor K appearing in equation (|34|), This unknown factor K is fixed by requiring that the total 
mass of the system, equation (|||), satisfy M = 1. Next we integrate the field equations (Q)-® for 
p, 7, and u>, and then equation (f20|) for a. We then iterate this procedure until some convergence 
criterion is met. 

Using a combination of integral and finite differencing techniques we solve the equations for 
the matter and gravitational fields. These equations are solved on a discrete grid in p and r on the 
computational domain < r < 00 and < p < 1. Unlike the Newtonian case, we cannot restrict 
the computational domain to the matter interior because the matter-independent effective source 
terms R p , i? 7 and i? w are nonzero in the vacuum exterior. Consequently we divide the radial grid 
into an interior and exterior domain. Each domain is covered by a geometrically spaced grid in r, 
with the grids joined smoothly at the cluster surface. We use an angular grid that is uniformly 
spaced in p. The interior radial grid is arranged to yield sufficient resolution for the core of the 
cluster, while the outer grid extends to some sufficiently large radius, typically 3 times the radius 
of the matter surface. High resolution of the core is crucial for obtaining numerical accuracy in 
highly centrally condensed relativistic disks. 

The three elliptic field equations (Q) - @ are solved by an integral Green's function approach 
following Komatsu et al. (1989) and pook et al. (199"2) . Again we make a matter-dependent, 



matter-independent split of the Green's functions depending on whether we are integrating the 
matter-dependent source terms t p , t 7 , and or the matter-independent source terms R p , i? 7 and 
R^. Thus we write the solution of equation (|7j) as 

p(r,p) = W p (r,p) + G p (r,p) (51) 

where we have defined 

00 POO fl 

W p (r,p) =-Y, e" 7/2 / dr' \ dp'r' 2 fl(r,r')P 2n (p)P 2n (p')R p (r' , p>), (52) 

n=0 



JO 



and 



G P (r,p) =~Y, ^ /2 P2n(p) [ j " I dr'r'f 2n (r,r%(r',0). (53) 

71=0 







Here we have used the 5-function in p to carry out the p! integration in G p . Similarly to solve 
equation (||) we write 

rsin$7(r, p) = W 1 (r,p) + G 7 (r, p) (54) 

with 

00 roo rl 1 

W 7 (r,p) = — J2e~^ 2 dr' dp'r' 2 ft-iM- sin(2n - 1)0 

7T ^ Jo Jo 2n - 1 
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and 



G y (r,n) = - E e 



xsin(2n- l)6>%(r',//), 
1 ^ _ 7/2 sin(2n-l)g(-l)" 



n=l 



(2n - 1) 

Finally, to solve equation (||) we write 

r sin 6 u(r,(i) = W w (r,p) + G u (r,p,) 

with 



dr'fL-i(r,r%(r',0). 



(55) 
(56) 

(57) 



,(2p~7)/2 



("OO /-l ^ 

Jo *7o ^^'^^' ^nVn-l) 



n=l 



and 



G w (r,/i) = -5 E 



,(2p-7)/2 



(-lf(2n-l)!! 



n=l 



Here 



and 



2n(2ra-l) 2«(n-l)! Jo 

fVr r'l = / ( r '/ r )"' for r '/ r ^ 1 
I (r/r'r, for r'/r > 1 



drVV£»-ifcr%(r',0). 



2 f (l/r)(r'/r)«, for r'/r < 1, 

Jn[r,r ' \ (l/r')(r/r') n , for r'/r > 1 



(58) 



(59) 



(60) 



(61) 



Among the advantages of this Green's function approach for solving the elliptic field equations 
is that the asymptotic conditions on p, 7 and lo are imposed automatically. That is, p ~ 0(l/r), 
7 ~ C(l/r 2 ), and uj ~ 0(l/r 3 ) for large r. To improve the accuracy of the angular integrations, 



we use the identities in equations (34) - (38) of Cook et al. (1992) 



It should be noted that in the integration of equation (^), the calculation of t u {r, 0) at r = 
requires special care. It can be seen from equations (|33|) , (|39|) and (f4G|) that tos(r = 0) = 0, while 
regularity conditions near the axis imply that £03 ~ r near r = 0. Thus, at r = 0, the calculation 
of the quantity t^^/r that appears in equation (|l9|) is done by computing the quantity dtos / dr\ r= o 
analytically. We do this by taking a derivative of equation (|33|) directly. For any r this gives 



, . , d 
dt 03 /dr = — 



pt (r) 



dp^p 4 



dp f f GR (E, J z 



(62) 



vlir) $ dl(p*, r ) 

dp v p v 

pt(r) or 



2 ? dj) if} Ay Ay dp (^*) 

+p+(r)I(pt-{r),r) — ±- p_(r)I(p_(r),r)- 



dr 



dr 



(63) 
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Here p±(r) and p r m ax{ r iP ) are defined in equations ( |40| ) and (^9|) and we have denned I(p^,r) to 
be 

r ) = f Pmaxir ' P ] dp f f GR (E, J z ). (64) 



o 



At r = it can be shown that the first term in equation (|63| ) vanishes by symmetry. The other 
two terms involve the integral I(p^, r) evaluated at p^ = p±(0) and r = 0. From equation ( |64|) and 
the form of the distribution function fen (equation fl34])) we see that this integrand is divergent 
and the limits collapse to zero at p& = p±(0). However, near p^ = p±(0) we can expand the 
integrand (and upper limit) to get an estimate of the value of the integral for use in the calculation 
of dto3/dr\ r= Q. This yields 



dt 03 /dr\ r=0 = 7r 



^a(r){a 2 {r) - 1) a(r) d r (r)e~^/ 2 



(65) 

r=0 



Here a(r) and d(r) are functions defined by equations (41) and ( |4^ ) 



Equation (^) for a is solved by integrating the linear ODE from the pole (fj, = 1) to the 
equator with the initial condition that 

a = at fi = l, (66) 

which arises from the requirement of local flatness on the coordinate axis. The derivatives of p, 7 
and to appearing in the matter-independent source terms R p , Ry, and the right-hand side of 
equation (^0|) are evaluated by finite differencing. 



Since we can evaluate the integrands in equations (pl|) - (|33|) at any values of p r and p®, we 



carry out the quadrature over each variable by Romberg integration ( Press et al. 1992 ) 



For a convergent solution, we require the maximum fractional change in all four metric 
functions to be less than 1 x 10~ 3 on successive iterations. This typically takes about 20 iterations 
when we start with a Newtonian initial guess (Section ^). A significant number of iterations can 
be saved during a sequence of calculations in which the heating parameter b is varied if we use the 
previous solution as the initial guess for the next value of b. 

For disks in the relativistic region, 600 radial zones, 20 angular zones, and 20 Legendre 
polynomials are adequate. 



6. Initial Data 

As an initial guess for the metric potentials p(r,fi), 7(7", (j), u>(r,p), cr(r,p) we used their 
Newtonian limits, given in terms of the Newtonian potential $at, by 



p(r,p) -> 2$ N (r,n), 
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7(r, fj) 
a(r, n) 



0, 



0. 



(67) 



Mihalas (1969) gives an expression for the Newtonian potential of an oblate homogeneous spheroid. 



By taking the limit as the eccentricity e — * 1, we obtain the Newtonian potential 3>jv(r, (j) for all 
points both on and off the disk. Outside the matter, we find 



*tf(r,/x) 



3M 

*2Rrn 



D 



2(1- X 2)l/2 



(1 - 3// 2 )(£(l - x 2 ) 1/2 - x) + (1 - ^ 2 )s 



where g s = r~ jR m , x = yfp, and i? = arcsinx with p defined by the quadratic equation 



Inside the matter, we have 



gi (l-^V-(g s + l)p + l = 0. 



3irM , q„ 
Mr,M) = -^(l-|). 



(68) 

(69) 
(70) 



7. Numerical Results 

Our region of investigation spans values of the heating parameter in the range .05 < b < .95 
for values of R m /M ranging from R m /M = 455 down to R m /M = 6.26. Figure |l| (a) shows the 
convergence of some representative runs in the Newtonian region. As we can see, as b increases, 
R m /M of the solution decreases. This is also evident for the more relativistic cases shown in 
Figure |l| (b) . 

To better understand our solutions it is helpful to look at how the central redshift, z c , varies 
with the heating parameter, b. For the Newtonian Kalnajs disk given in Appendix A, we can 
calculate the central redshift, z c analytically from the Newtonian gravitational potential as follows: 

z c = [-(teCO.O)]" 1 / 2 -! (71) 



e 7(0,0)+p(0,0) 



-1/2 



Now use equations (|67j) to get 



^Newtonian = _^( 0) Q) = ¥£L. (72) 

AR 



This quantity is independent of the heating parameter b and thus is a useful diagnostic. Although 
analytic equilibrium solutions exist in the Newtonian case even in the unstable region (6 > .816) 
(see Appendix A), it is not obvious that our iterative method will converge to them. Instabilities 
may be signaled by failure to converge, or by causing convergence to a solution that is far from the 
true Kalnajs solution. Thus by plotting the central redshift versus the heating parameter we can 
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Fig. 1. — Convergent runs for some representative values of R m /M. (a) Runs in the Newtonian 
region, (b) Runs in the more relativistic region. 

hope to see where the solution we find differs from the true Newtonian solution. As can be seen in 
Figure || (a), for a Newtonian cluster with R m /M = 455, our code finds the Kalnajs equilibrium 
solutions for b < .6. Since the stability region is given by < b < .816, the change in the central 
redshift does not serve as a sharp indicator of the transition to the unstable region. 

Similar behavior can be observed in our relativistic clusters. Figure |^ (b) shows an example 
for R m /M = 8.84. Regions in which z c is not constant may represent true general relativistic 
equilibria, but most likely these clusters are near or inside the general relativistic unstable region. 
To fully pin down this region requires a dynamical stability analysis, which can probably only be 
done by numerical evolution. 

The most relativistic solution we have found has an R m /M of 6.97. By comparing this disk 
to the Kerr geometry we can get an idea of how relativistic this model is. Using the proper 
circumference of the disk we find that our most relativistic cluster corresponds to an R/M for 
Kerr of 8.09. Table | gives some properties of the representative models discussed here. 

In Figure || we show contour plots of the surface density for two examples of our equilibrium 
disks. Frame (a) shows our most Newtonian cluster (R m /M = 455) with b = .6. In Newtonian 
theory these disks have homogeneous volume density. Frame (b) shows the same plot for the 
analytic version of the same cluster. In contrast, Frames (c) and (d) show our most relativistic 
cluster (Rm/M = 6.97) with b = .1 and its Newtonian counterpart respectively. The relativistic 
disks are in general more centrally condensed than their Newtonian cousins. 
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Fig. 2. — Central redshift {z c ) versus heating parameter b for our most Newtonian cluster, 
(a) R m /M = 455 and a representative relativistic cluster, (b) R m /M = 8.84 . 



In previous work ( Shapiro Sz Teukolsky 1985,1992b| ) on polytropic axisymmetric clusters we 



found that a maximum in the binding energy along an equilibrium sequence marked the onset of a 
dynamical instability. Since, for the Newtonian Kalnajs disks here the potential, (equation ( |A7[ )) 
and surface matter density, (equation (|A1])) are not functions of the heating parameter, the 
binding energy is also not a function of b. Thus it cannot be used to analyze stability. However, 
in the general relativistic case we expect this degeneracy to be lifted. We do indeed find that the 
binding energy does depend on b. However, because of numerical errors, we are not able to reliably 
find the turning point in the binding energy that would signal the onset of instability. Again, a 
firm limit on the stability region must wait for a full dynamical evolution code. 
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A. Newtonian Kalnajs Disks 

A Newtonian Kalnajs disk can be obtained by "flattening" a homogeneous oblate spheroid, 
i.e., by letting the eccentricity e — > 1 ( |Binney Sz Tremaine 1987| ). This flattening yields a disk with 
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Fig. 3. — Surface density contour plots: (a) Our most Newtonian cluster with R m /M = 455. and 
(b) its analytic analogue, (c) Our most relativistic cluster with R m /M = 6.97. (d) The analytic 
Newtonian version of this cluster {R m /M = 6.97). 

a surface density given by 

Z(R) = Z C (1-R 2 /R 2 J, (Al) 

where S c is the central density, R m is the radius of the disk (matter radius) and R is the radius in 
the disk plane defined by 

R = rsm9. (A2) 

When all the particles are in circular orbits with uniform angular velocity fi C i rc , is 
related to the central density by the same relation as for Maclaurin spheroids. 
I . 1 

Kalnajs (1972) developed a family of equilibrium disks, all of which have the above surface 



density. Each member of this family is characterized by the parameter b defined in equation (35) 
These disks have velocity dispersions that are governed by the tunable parameter b: 



(vi) 2 - U\ = (yR) 2 = V(l _ - Rt/Rl). (A3) 



3 

From the above equation we can see that when b ~ 1 {QR m ~ V) these disks are "cold", i.e. 
have very little thermal motion. This corresponds physically to a system in which all particles 
move on nearly circular orbits. On the other hand, when 6^1 {QR m <C V) we have "hot" 
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systems in which most of the support against self-gravity comes from the random motions given 
in equation QA3| ) ( [Binney fc Tremaine 1987| ). 



A distribution functions that generates this family is 



In (En, J z n) 



K N 




2{(bVJ zN /R m ) - E N } - V 2 (l + b 2 ) 



-1/2 



[-]>0 

[•••]<o 



(A4) 



Here Kn is a constant that can be determined by fixing the total mass of the system, En is the 
Newtonian energy, and J z n is the Newtonian angular momentum about the symmetry (z) axis. 
These quantities, with m = 1, are given by 



En = — ^ 1 ~ h 4>n(R) 



and 



JzN = Rp^- 

The Newtonian potential, <f>^(R), in the disk plane is 

N (R) = <S> N (r,O) = -V 2 (l 



2R 2 ' 



(A5) 
(A6) 

(A7) 



Using the linearized collisionless Boltzmann equation (linearized Vlasov equation) Kalnaj 



(1972) analyzed the stability of this family of disks. He found that restricting the heating 
parameter, b, to be less than 0.816 produced disk systems that are stable against all axisymmetric 
disturbances, i.e. disks with < b < .816 are stable against ring formation. 
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Table 1. Properties of some representative solutions. 



Rm/M 


b 


Zc 


„max a 
9tt 




J c 


455 


0.6 


0.0051 


-0.9897 


-0.0032 


4.2304 


6.97 


0.1 


.9556 


-0.2614 


d 


0.0693 



& 9tt ax is the maximum value of gu- 
b Ei, is in units of M . 
c J is in units of M 2 

d For the relativistic disk here our value of E& is 
unreliable because of numerical errors. 



